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In a recent review it is said that free-surface flows “represent some of the difficult 
remaining challenges in computational fluid dynamics” [Tsai and Yue, 1996]. There has 
been progress with the development of new approaches to treating interfaces, such as the 
level-set method and the improvement of older methods such as the VOF method [Kothe 
et al., 1996]. 

A common theme of many of the new developments has been the regularization of 
discontinuities at the interface. One example of this approach is the continuum surface 
force (CSF) formulation for surface tension, which replaces the surface stress given by 
Laplace’s equation by an equivalent volume force [Brackbill et al . , 1992]. 

Here, we describe how CSF formulation might be made more useful. Specifically, 
we consider a derivation of the CSF equations from a minimization of surface energy 
as outlined by Jacqmin [1996]. This reformulation suggests that if one eliminates the 
computation of curvature in terms of a unit normal vector, parasitic currents may be 
eliminated. For this reformulation to work, it is necessary that transition region thickness 
be controlled. Various means for this, in addition to the one discussed by Jacqmin [1996], 
are discussed. 

Derivation of the CSF Equations from an Energy Functional 

Beginning with the energy attributed by Jacqmin [1996] to van der Waals, we show 
that minimization of the surface energy is equivalent to Plateau’s problem [Logan] , in which 
surface area is minimized. Further, following Jacqmin’s virtual work argument, we derive 
a CSF-like formulation. The difference between this and the standard CSF formulation 
illuminates the origin of parasitic currents and suggests that these are best eliminated by 
controlling the interface thickness. 

Consider an incompressible fluid with density pi in region fl separated by a closed 
surface dfl from a second incompressible fluid with density p 2 \ 


f pi , xCfi 

\ P2 , X £ H . 


( 1 ) 


On dQ, in the usual surface tension formulation given by Laplace’s equation, a surface 
stress boundary condition is imposed, 


(Pi -P2+ CTK)h = 0 , 


( 2 ) 


where P 12 is the pressure in fluids 1 and 2, a is the surface tension coefficient, and k is 
the local surface curvature. 

In the CSF model, the density, p, is replaced by a regularized density that is calculated 
by convolving p(x) with an interpolation function, S’(x / — x;/i). S has bounded support 
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h, is normalized, and decreases monotonically with increasing |x — x'|. The regularized 
density, 

p(x) = p J P(x')5(x - x')d 3 x' , (3) 

is differentiable, 

Vp(x) = ^ J p(x')ViS(x - x>V . (4) 

Consider the term in the surface-energy integral, given by Jacqmin [1996]. 

e = J d 3 x ^aVp(x) • Vp(x) + /3/(p) , (5) 

I 

in which the first term is proportional to composition gradients. As in Brackbill et al. 
[1992], Vp is given by an integral over the interface, 

Vp(x) = W / n(x»)5(x - x s )<M , (6) 

h Jdci 

where [p] — p\ — p 2 - Substituting this expression into Eq. (5) and setting /? = 0 gives 

I = TT f [ n(x s ) • n(x' s )T(x s ,x';/i)dAdA' , (7) 

« Van Van' 

where T is a “mass-matrix” defined by 

T(x s , x' s ; h) = ^ J d 3 xS(x - x' )S(x - x.) . (8) 

Clearly, when h is small compared with the local radius of curvature, T — S(x 3 — x' s ) , and 
the ratio of e to the area of the interface is a constant, 

A = e/[p] 2 /h 2 . (9) 

Thus, minimizing the energy integral (with /3 = 0) is equivalent to minimizing the surface 
area, a classic variational problem called Plateau’s problem [Logan, 1987]. 

Following Jacqmin, we include the surface energy contribution in the total energy of 
the fluid, 

E = J±pu 2 + ±a(Vp) 2 dV . (10) 

For simplicity, we neglect the wall-surface energy. Energy is a constant of the motion, so 
that 

<u) 

The first term is the work done by forces acting on the fluid, 

|^ u2 = u p < 12 > 
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The second term can be integrated by parts to yield 

0 = J u • [F + a (V 2 p) Vp] dV , 

where we have substituted the continuity condition for an incompressible fluid, 

dp 


(13) 


dt 


= — u • Vp . 


(14) 


Conservation of energy must hold whatever the velocity. Thus the surface tension force is 
given by 

F = —a (V 2 p) Vp 


(15) 


or substituting n = Vp, 


By comparison, the CSF force is 


F = —a (V • n) n . 


Fcsf = <tk(x) 


Vp(x) 

M ’ 


where the curvature is given by 
Since n = -^n, then 


k — —V • n . 


V • n = ^V-n + 0 (h 

h 


R 


and Fcsf and the variational force are equal when 


a 


h 

\P\ 


-a . 


(16) 


(17) 


(18) 


(19) 


(Jacqmin, using the incompressible flow condition, 

V • u = 0 , 


( 20 ) 


integrates by parts the second term in Eq. (13) to derive the potential form of the surface 
tension, 

F = -pVp , (21) 

where p = — oV 2 p is the chemical potential.) 


Parasitic Currents 

Parasitic currents are a troubling artifact of the CSF formulation. These slowly 
growing vortical flows in the transition region are driven by a vorticity source term equal 
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to V x F /p. Jacqmin’s form, Eq. (21), gives no contribution to this term. Equation (15) 
does give a contribution, 

V x F/p = --V(V 2 p) x Vp . (22) 

P 

T his contribution will be nonzero until the “curvature,” V 2 p, of the interface is equal to a 
constant. Since this minimizes the surface energy, it would seem to be the correct physical 
result. The CSF formulation gives a similar contribution, 

V x F/p = -j— ^ Vk x Vp . (23) 

\P\ 

However, the calculation of the unit normal from the gradient of the density requires some 
care to avoid singularity, which appears to cause relatively large alterations of the curvature 
with small displacements of the interface. These, we believe, are the most probable source 
of the parasitic currents. 


Comments on Controlling Transition Region Thickness 

The essential difference between the CSF formulation, which was derived by the 
construction of a surface delta function that reproduces Laplace’s equation in the limit 
of vanishing h , and the variational form is in the expression for the curvature. In the 
variational form, V • n replaces k = -V • n. This substitution is valid if the transition 
region remains narrow. 

Particle-in-cell methods and front-tracking methods in general provide solutions which 
preserve narrow transition regions. Further, as Kothe et al. note [2], the VOF technique 
eliminates numerical diffusion normal to the interface. (An example of a VOF calculation 
in three dimensions is shown in Fig. 1.) Thus, it is reasonable to expect that a reformulation 
of the CSF model which eliminates the calculation of the unit normal can be successfully 
implemented in a numerical code. 

One can also consider techniques that actively control interface thickness. For 
example, in Jacqmin’s formulation, the second term in the surface energy integral, Eq. (5), 
controls the thickness of the interface. Jacqmin notes that its effect is to segregate two 
immiscible fluids. By appropriate choices of /(p), this term introduces a nonlinear diffusion 
that tends to narrow the transition region. 

The use of this term in the energy integral to steepen the interface appeals to a 
physical cause. However, interface steepening can be introduced in a formulation that is 
purely numerical. Yabe [Yabe and Xiao, 1995], for example, discusses the convection of 
the density, Eq. (1), described in one dimension by 


dp dp 
dt Ul dxi 


= 0 . 


(24) 


The same equation describes the advection of any function of p, 


dFjp ) dFjp) 

dt Ut dx t 


= 0 . 


(25) 
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Specifically, Yabe and Xiao consider the tangent function, 


F(p) = tan [0.997T (p/{p) — 0.5)] , 


( 26 ) 


where (p) = l/2(pi + p 2 )- The factor 0.99 avoids infinite values when p = pi or p 2 . It is 
necessary, of course, that p\ < p < p 2 . 

The technique maintains the sharpness of the transition region in the following way. 
The numerical solution of Eq. (25) will result in some numerical diffusion with diffusivity, 
a = a (Ax n , A t, m ). Thus, the associated equation can be written (in one dimension) 


dF_ dF _ d 2 F 
dt Ul dxi ° dxf 


(27) 


Substituting this equation into the original advection equation (24), with the transforma- 
tion given by Eq. (26), one finds that solutions of the ordinary differential equation, 


are stationary solutions of Eq. (24). The solution to this equation, 

p = - tan" 1 (e**) , (29) 

7T 

is a continuous function with a steep transition between two nearly constant states. 
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Figure Captions 

A shallow-pool drop splash calculation is shown in Figs. 1-4. The mesh spans [0.0, 8.0] 
(cm) in the x and y directions and [-8. 0,8.0] (cm) in the z direction with 16 x 16 x 32 cells. 
A spherical drop with a 4.0 cm radius is centered at (0.0, 0.0, 3.0) above a shallow pool in 
the xy plane that occupies 2 = [—8.0, —5.5]. The drop and pool are the same fluid, with 
10 times the density of the background fluid density. Surface tension is neglected and the 
fluids are assumed to be inviscid. Gravity is (0.0,0.0,-980.6) cm/s 2 . Interfaces are tracked 
with our new 3-D piecewise planar volume tracking algorithm, and the incompressible flow 
is modeled with an approximate projection technique. 

The drop size is large compared to the shallow pool (having a volume 25% of the pool 
volume), so the drop impact on the pool generates a sizable splash and rebound jet, as 
seen in Figs. 3 and 4. Plotted are isosurfaces of the 1/2 volume fraction level at t = 0 (1), 
t = 0.12 (2), t = 0.20 (3), and t = 0.40 (4). 
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